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Abstract 

We describe a Monte Carlo procedure for the simulation of dy- 
namically triangulate random surfaces with a boundary (topol- 
ogy of a disk). The algorithm keeps the total number of triangles 
fixed, while the length of the boundary is allowed to fluctuate. 
The algorithm works in the presence of matter fields. 

We here present results for the pure gravity case. The algo- 
rithm reproduces the theoretical expectations. 



1 Introduction 



The quantization of the 2-dimensional Einstein-Hilbert theory of gravity 
presents conceptual and technical difficulties In the Euclidean path in- 
tegral formulation, one has to integrate over all metrics (modulo diffeomor- 
phisms) of a Riemannian manifold with fixed topology. One way of making 
sense out of this ill defined integral is to approximate the manifolds by trian- 
gulations. The integration over all metrics is then replaced by a summation 
over all triangulations. 

The approach to quantum gravity via Dynamically Triangulate Random 
Surfaces (DTRS) opens the way for the use of Monte Carlo methods which 
were proven to be very useful, especially if it comes to the inclusion of matter 
fields [|]. 

In this letter, we present a Monte Carlo procedure for the simulation of 
DTRS that have a boundary. 

The importance of the boundary comes from the fact that the estimation 
of the probability distribution of its length corresponds to a measurement 
of the wave function of the 2D universe. Its theoretical treatment turned 
out to be problematic. With only the Hilbert-Einstein term in the action, 
2-D quantum gravity is non-renormalizable by naive power counting and is 
known to suffer from ultraviolet divergences. Furthermore, when conformal 
field theory is coupled to 2-D gravity, a naive counting yields "—1" degrees 
of freedom 0. This is believed to reflect a non-normalizable universe wave 
function solution of the Wheeler-de Witt equations. Matrix model results 
indicate too that the wave function might be non-normalizable at small areas. 
These considerations led us to a numerical study of open random surfaces 
(surfaces with a boundary). The hope is that adding extra boundary terms 
(or additional matter flelds) to the action might lead to a normalizable wave 
function. 

It is therefore useful to have a Monte Carlo algorithm for the simulation 
of open random surfaces, not only to understand the nature of a possible 
normalizable wave function but also as a tool to study the intrinsic geometry 
of open surfaces. 

A Monte Carlo procedure for the simulation of closed surfaces with a 
flxed number of triangles, using the (2,2) Alexander move (to be called "bulk 
flip" later) was studied by Kazakov et al. A generalization of this type 
of algorithm to the case of a surface with a boundary with fluctuating length 
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turns out to be nontrivial. The main problem is to ensure ergodicity and de- 
tailed balance (i.e. that each triangulation contributes with equal probability 
to the partition function). 

This letter is organized as follows. The model is introduced in section ^ 
Section § is devoted to the description of the algorithm. In section |^ we 
present our results and compare them with the theoretical values. 

In the present paper we limit ourselves to the pure gravity case where 
most of our results can be obtained also analytically. The inclusion of the 
matter fields (possibly with central charge c > 1) where no analytical results 
are available is straightforward algorithmically, but its theoretical diagnostics 
and interpretation is beyond the framework of the present letter. 



2 The model 

The Euclidean Einstein-Hilbert action is 

SE[9] = -^f d'x^\{-R + 2A). (1) 

Here is Newton's constant, R is the intrinsic curvature, and A is the 
cosmo logical constant, in units oi h = c = 1. Misa Riemannian space time 
manifold. 

In the continuum quantized version (for the sake of completeness we in- 
clude matter fields X here) one considers the partition function 

The discrete analog of eq. (H) without matter fields is 

^ = E ^ exp (aoA^o + aiA^i + asA^s) ■ (3) 

The sum is over all onentec? triangulations r of the surface. No,Ni,N2 are 
the number of vertices, links and triangles of r, respectively, and s(r) is 
the order of the symmetry group of the triangulation r. ao,ai,a2 are free 
parameters. The factor can be viewed as what remains from the ^^^^^.^^ 
factor in eq. (H). 
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For open surfaces the action in eq. (0) picks up two additional terms (as 
in the Liouville theory), 

^^^ + T7r?^f ds\/h{-k + 4:E), (4) 
IottGat JdM 

where h is the induced boundary metric, k is the extrinsic curvature on the 
boundary, and S is the boundary cosmo logical constant. 

The Einstein-Hilbert action is a topological invariant, since the 2-D ori- 
ented manifolds obey the Gauss-Bonnet theorem (with the correction arising 
from the boundary term): 

f d^x\f\g\R+(( \[}idsk)=x- (5) 

A-K Jm ^ JdM 

Here x = 2 — 2/i — 6 is the Euler characteristic, h is the number of handles, 
and b is the number of boundaries. Nevertheless some non-trivial quantum 
theory exists, mainly due to conformal anomalies f^. 

From the topological relation x = A^o — A^^i + A^2, and the relation for a 
connected triangular lattice 2Ni = 3N2 + L, where L is the boundary length, 
the discrete action in eq. (|^) can be written as 

aoNo + aiNi + a2N2 = -XN2 - + 7^ , (6) 

with 

-i- + — + a2) 
^2 2 ^ 

ao + ai 
2 

ao. (7) 

The action depends on three independent parameters: A and ^, that we 
identify with the cosmological constants A and H, and 7, that we identify 
with I/Gat. Actually as we fix x, since we identify 7 with a constant term, 
we expect the partition function to look like 

Z = Y,^ e-"^ e-«^ = J2 Z{L, N) e-^^e'^^ . (8) 

{r} ^^''~) N,L 



A = 

e = 

7 = 
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Here and in the following we have identified N = N2. 

In our simulations to be specified below one estimates Z{N, L) which is 
the combinatorial weight of a surface with N triangles and L external edges. 



We consider open connected triangulate surfaces without 1- and 2-loopsE 
In the language of the dual lattice, where a triangle becomes a point, an 
edge a segment and a vertex a polygon, every surface 5* is in one-to-one 
correspondence with a connected planar diagram of the (p^ theory without 
tadpoles and self-energy. 

One of our aims is to find the wave function ip{L)^ which is the probabil- 
ity amplitude to have a boundary length L, for a fixed number of triangles 
N . Our Monte Carlo procedure samples the configuration space with no a 
priori weight (each configuration occurs with the same probability). We can 
therefore determine the probability that a configuration has boundary length 
L by just making a simple histogram, cf. section ^. 

As a by-product of the simulation one gets knowledge about the intrinsic 
geometry of the surface. Using the simulation results, we are able to picture 
and characterize the typical surfaces. 

3 Algorithm 

We consider a connected, open, triangulate lattice, with the topology of a 
disk. By definition, all edges have equal length. We have to distinguish three 
types of triangles: 

• type-0 triangles have none of their three edges on the boundary of the 
surface (bulk triangles), 

• type-1 triangles share exactly one edge with the boundary, 

• type-2 triangles have two edges in common with the boundary. 

Our Monte Carlo algorithm generates a start configuration that consists of 
N triangles and that has a certain boundary length Lq. Then a number of 
sweeps is done. Each sweep consists of a sequence of elementary fiips done 
on the configuration: 

^ Given a triangulate surface S, a p-loop is a set of p distinct edges in S which form a 
loop on S 
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bulk flips change the internal geometry but maintain a fixed boundary 
length. 



• boundary flips change the boundary length L. There are two kinds of 
boundary fiips: 

— type-1 boundary flips are performed on type-1 triangles and in- 
crease the boundary length by 2. 

— type-2 boundary flips act on type-2 triangles and decrease the 
boundary length by 2. 

Because of the relation E = ^^2^ for a connected triangular lattice, where 
E is the (integer) number of edges in the configuration, we deduce that for 
fixed iV, the boundary length L may be incremented by even values only. 

Let us characterize in detail each of the three fiip operations: 

• Bulk flips (see fig. |I]): 

These are the basic moves also used in the simulation of closed DTRS 
i, see fig. |. 

• Type-1 boundary flips (see fig. |^): 

One adds a vertex to the configuration and connects it to the two border 
sites of the type-1 triangle. The connecting links are called "virtual" 
edges. One then performs a bulk flip on the former border link of the 
type-1 triangle. Now one selects one of the two virtual boundary links 
(each of them with probability one half) and removes it. The remaining 
virtual link is then made "real". The so defined operation obviously 
increases L by 2. Note that the type-1 fiip is allowed only when the 
resulting configuration is connected. 

• Type-2 boundary flips (see fig. |^): 

To the given type-2 triangle one adds two virtual links as depicted in 
fig. ^ Then the type-2 triangle is removed and, with uniform proba- 
bility, one of the two virtual links is selected and removed. Finally the 
remaining virtual link is made real. The type-2 boundary flip decreases 
L by 2. 

In order to be able to prove later the detailed balance condition (9), the 3 
types of fiips have to be applied according to the following prescription. 
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During the Monte Carlo sweeps one selects again and again with uniform 
probability one of the triangles as a candidate for updating. 

Once a triangle is selected, one decides what kind of flip will be attempted 
on it. 

Strange as it may look, the constraints of the detailed balance proof (see 
below) require that the kind of the update (flip) is to be decided a priori, 
independent of the actual type of the triangle on which it is going to be 
attempted. 

More precisely, one decides with probability < p < 1 to perform a 
boundary flip, and with probability 1 — p to perform a bulk flip {p is a free 
parameter which we chose in our runs to be 0.5).i 

If the selected update kind is "boundary" and the triangle is type-0 
("bulk") then the present updating attempt is aborted and one proceeds 
to the selection of a new candidate triangle. 

If the update kind is "boundary" and the triangle is type-1 or type-2 then 
one applies the corresponding boundary flip. 

If the update kind is "bulk" one selects with uniform probability one of 
the three links (edges) of the triangle. If the link belongs to the boundary, 
the updating attempt is aborted and one proceeds to the selection of a new 
triangle. Otherwise the selected link is bulk flipped. 

In Monte-Carlo simulations one must ensure that the stochastic process is 
ergodic, i.e., that all configurations can be reached starting from an arbitrary 
initial configuration. 

The ergodicity for fixed and L comes from the ergodicity of the Alexan- 
der (2,2) move (the bulk flip). Therefore all configurations with fixed {N,L) 
can be reached by a finite number of steps. With boundary flips one can get, 
for a fixed A^, from any L to any L provided that L and L have the same 
parity and that 3 < L < A^ + 2. These two restrictions are specific for the 
model that we consider. Thus the combined flips provide an ergodic process 
for a fixed number of triangles and a varying boundary length. 

In order to ensure that each configuration is generated with equal prob- 
ability we demand that our Monte Carlo procedure satisfies the detailed 
balance condition. In the absence of matter fields this means that the prob- 

^Another flip operation, also consistent with stationarity, is performing bulk and bound- 
ary flips in fixed order, e.g., alternating 
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ability P to move from one configuration to another obeys 

P(r ^ f ) = P{f t) . (9) 

It is known that eq. (|^) is fulfilled for bulk flips We therefore consider 
configurations f that differ from r by a boundary flip. The probability that 
a certain triangle is selected to be boundary flipped is j^. (Of course, as 
mentioned above, the triangle will be actually flipped, only if is not a type-0 
triangle.) There are four flip possibilities. In each one a couple of triangles 
(*,*) from T (where (*,*) expresses the types of the triangles) change their 
type: 

(1+) a type-1 boundary flip changes a {type- 1, type-0) combination into a 
{type- 2, type- 1) combination. This operation is applied with probability 
2^. The factor ^ comes from choosing the neighbour type-0 triangle 
out of two possible neighbours of the type-1 triangle. 

(1— ) a type-2 boundary flip (also selected with probability moves a {type- 
l,type-2) combination to a {type-0, type-1) combination. This is the 
move depicted in fig. ^. It is the "inverse" operation to (1+). 

(2+) a type-1 boundary flip (selected with probability j^) changes a {type- 
1, type-1) into a {type-2 , type-2) combination. This is the operation 
depicted in fig. 

(2—) The "inverse" operation of (2+) is a type-2 boundary flip taking a {type- 
1, type-2) combination into a {type- 1 , type- 1) combination. It is also 
selected with probability j^. 

The verification of the detailed balance condition is in the statement that 
the (+) and (-) operations are inverse to each other and are selected with 
the same probability. 

Note that our algorithm can be easily adapted to the inclusion of matter 
fields (which live on the triangles). 

For the actual simulation we used two geometrically very different initial 
configurations. This enabled us to check the ergodicity of our algorithm and 
also gave us the possibility to check for critical slowing down. The two initial 
configurations were: 
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1. A polyhedron consisting of triangles, see fig. ^. All triangles are 
type-1 and share a single vertex. 

2. A configuration that contains 20 type-1 triangles and A^ — 20 type-0 
triangles. For large A^ this initial configuration resembles a sphere with 
a hole in it. 

This second initial configuration is constructed by taking a configuration 
of the first type with A^ = 20 and growing it by inserting vertices in the 
middle of the triangles 0, see fig. ^ Each vertex insertion increases A^ by 
two. After each insertion of a vertex, we performed a sweep. 

4 Results 

Using the algorithm described above we looked at various quantities. For sev- 
eral different values of A^ we obtained estimates for the mean of the boundary 
length {{L)), and the variance of the boundary length ((L^) — {L)'^)- We also 
studied the expectation values of the number of type-0, 1,2 triangles and of 
the number of triangles Nf that do not have any vertex on the boundary. 

For thermalization we typically performed 15,000 sweeps. For measure- 
ment we typically performed 150,000 sweeps per run (measuring after each 
sweep). 

(L) turns out to be a linear function of A^ with high precision. In table |I| 
we present our Monte Carlo results for A^ ranging from 50 to 6400. We give 
always two estimates for (L), one obtained from a run with initial boundary 
length Lq = 20 and one from a run with Lq = N. The results are always 
nicely consistent within la error bars. Fitting the data with the law (L) = 
A + BN we Gold A = 1.61(3), B = 0.764697(34), xV^of = 1.5 for the 
data obtained from runs with Lq = 20, and A = 1.62(3), B = 0.764643(34), 
X^/dof = 2.4 for the data obtained from runs with Lq = N. 

The constant B can be obtained analytically: The result for our lattice 
model is B = if = 0.764706 0. Our fit result from the Lq = 20 runs are 
nicely consistent with this exact result while the fit to the Lq = N result is 
consistent within 2a. 

The linear dependence of (L) on A^ indicates a fractal dimension of the 
typical equilibrium configuration and the fact that it is composed by loosely 
connected subgraphs. Further evidence that this is indeed the situation can 
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be seen in table ^ that quotes the expectation values of the number of tri- 
angles not touching the boundary (Nf). This data indicates that a typical 
configuration is tree-shaped with most links belonging to the boundary. Most 
of the branches look like beads of minimal width (1 link across). 

The variance of L goes like N2 . This suggests that the locations of the 
further branching within a branch are random independent variables. In 
particular, this implies that there are no long range correlations between the 
various local geometric features of these "surfaces" . 

To obtain further information about the global structure of equilibrium 
configurations we made some "snapshots" of small configurations. The 
"portrait" of such a typical surface (represented on the dual lattice) is shown 
in fig. 1^. This "snapshot" substantiates the picture that the typical surfaces 
are tree or root like with branching occurring at random walk steps on the 
configuration. 

Fig. shows an example for the statistical distribution of the border 
length L for = 400 (the wave function iIj{L)). The diamonds show the 
exact result , the points with error bars show estimates produced with our 
algorithm. The error bars were obtained as follows: The total of 100,000 
measurements of L was produced in 10 independent runs, each of size 10,000 
sweeps. The sweeps consisted of 10,000 flips of each type. A histogram was 
determined for each subsample separately. The error bars were obtained as 
the la error from the 10 independent measurements of the histogram. 

In conclusion we have presented results from a simulation study of pure 2- 
D quantum gravity on open DTRS's. The numerical results agree with results 
obtained by analytical calculations . In addition our simulation algorithm 
also allowed us to study the geometry of the typical surface configurations. 

Our results encourage further advance by adding matter fields. 
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Table 1: Monte Carlo estimates for (L) 



Figure 1: The hulk flip 
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(a) (b) (c) 

Figure 2: How the type-1 boundary flip works: (a) A vertex is 
added to the configuration and connected to the two border sites 
of the type-1 triangle, (b) A bulk flip is performed, (c) One of 
the two virtual boundary edges is removed. The remaining edges 
are made real. 




Figure 3: How the type-2 boundary flip works: (a) Two virtual 
edges are added to the configuration, (b) The type-2 triangle is 
removed, (c) One of the virtual links is removed, the other is 
made real. 
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Figure 4: Polyhedron start configuration 




Figure 5: Growing the initial configuration 
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Table 2: Monte Carlo estimates for (Nf) 
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Figure 6: Snapshot of a configuration on the dual lattice 
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Figure 7: Monte Carlo (points with error bars) and exact results 
(diamonds) for the wave function ip{L) for N = 400 
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